{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Pulse Design Tool\n", "\n", "This notebook demonstrates how to use the **FreeGSNKE Pulse Design Tool** (FPDT). \n", "\n", "The FPDT combines the **evolutive equilibrium solver** introduced in previous notebooks with a **virtual plasma control system (PCS)** that governs the time evolution of the plasma. It enables **closed-loop** plasma scenario and control design studies.\n", "\n", "The virtual PCS generates voltage requests for the active poloidal field coils, which are then passed to the evolutive solver. At present, FPDT is primarily intended for use during the **flat-top phase** of a discharge. While it may function during ramp-up or ramp-down phases, its performance in those regimes is not guaranteed.\n", "\n", "The **controllers** implemented in the virtual PCS are inspired by the MAST-U PCS and allow control of the following quantities:\n", "- Plasma current.\n", "- Geometric shape parameters (derived from the equilibrium).\n", "- Vertical plasma position.\n", "- Coil activation times.\n", "- Coil current and voltage limits.\n", "\n", "For **more detailed information on these controllers and features**, please refer to the **FPDT paper on arXiv**.\n", "\n", "The FPDT can be used to:\n", "- Plan challenging new plasma scenarios.\n", "- Test new control schemes and parameter settings.\n", "- Serve as a training tool for plasma experimentalists.\n", "\n", "In the following sections, we demonstrate how to set up the FPDT for the **controlled evolution of a MAST-U–like plasma**.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate a starting equilibrium\n", "\n", "As in prior notebooks, we need a starting equilibrium. \n", "\n", "This equilibrium will be the initial condition from which the FPDT will evolve the plasma." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# package\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")\n", "\n", "\n", "# initialise equilibrium\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # Radial range\n", " Zmin=-2.2, Zmax=2.2, # Vertical range\n", " nx=65, # Number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=65, # Number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", ") \n", "\n", "# initialise profiles\n", "from freegsnke.jtor_update import Lao85\n", "profiles = Lao85(\n", " eq=eq, # equilibrium object\n", " Ip=7.59e5, # plasma current\n", " fvac=-0.52, # fvac = rB_{tor}\n", " alpha=[362685, 17696], # Lao profiles parameters\n", " beta=[0.103, 0.753], # Lao profiles parameters\n", " alpha_logic=True,\n", " beta_logic=True,\n", ")\n", "\n", "# set initial coil currents (passive structure currents start at zero here)\n", "coil_currents = np.array([2766.10782056, 193.12768191, 4019.86873871, 4857.42357466,\n", " -726.99282376, -1696.92521166, -77.7252396 , 266.38207493,\n", " -73.01483716, -4169.86874271, -4518.28147417, 5.2 ])\n", "for i, key in enumerate(tokamak.coil_names[0:12]):\n", " eq.tokamak.set_coil_current(coil_label=key, current_value=coil_currents[i])\n", "\n", "# initialise the static solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# solve for the equilbirium\n", "GSStaticSolver.solve(\n", " eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-8,\n", " )\n", "\n", "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=70)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Prepare inputs for the PCS class\n", "\n", "Here, we will set up the inputs (lists and dictionaries) for each of the (internal) controller classes within the virtual PCS (e.g. plasma, shape, virtual circuits). \n", "\n", "We start by defining which coils have which purpose." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# COIL NAMES AND SHAPE TARGET NAMES\n", "\n", "# active coils: must match the order they appear in FreeGSNKE machine description\n", "active_coils = ['Solenoid', 'PX', 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5', 'P6']\n", "\n", "# control coils: a subset of the active coils to be used for plasma and shape parameter control (i.e. we exclude vertical control coil P6 here)\n", "ctrl_coils = ['Solenoid', 'PX', 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n", "\n", "# vertical control coil\n", "vertical_coils=[\"P6\"]\n", "\n", "# names of the plasma targets (this will be clear later, equal to the number of solenoid coils)\n", "plasma_targets = ['plasma']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All of the following inputs are going to be specified as **waveforms**, i.e. they will be a dictionary of **times** and **vals**. \n", "\n", "For each array/list of **times** (in seconds), we must specify a corresponding list of list/arrays in **vals**. \n", "\n", "Most waveform quantities (e.g. `*_ref`, `*_ff`, `blend`) will use **linear** interpolation inside the PCS class internally, to enable querying at a given time $t$ within the simulation.\n", "\n", "Some other quantities (e.g. PID gains, virtual circuits, and other arrays), will use **previous-value (or step)** interpolation where the value defined at time `tmin` is used over the interval $ [t_\\mathrm{min},\\, t_\\mathrm{next})$, after which the value defined at `tnext` is applied.\n", "\n", "We will denote which quantity uses which type of interpolation with **linear** or **step**. \n", "\n", "The following two cells provide an example of how to set these two different waveforms types and what they look like when plotted. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# waveform dictionary structure (for linearly interpolated quantity)\n", "# e.g. this could be the desired radial position of the X-point over time\n", "waveform = {\n", " \"times\": np.array([0.0, 0.05, 0.15, 0.35, 0.45, 0.5]),\n", " \"vals\": np.array([0.57, 0.57, 0.5, 0.5, 0.57, 0.57]),\n", "}\n", "\n", "# plot what it looks like \n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "ax.plot(waveform['times'], waveform['vals'], color='k', linewidth=1, linestyle=\"--\", marker=\"x\", markersize=7, label=\"linearly interpolated waveform\")\n", "\n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Scalar quantity [units]\")\n", "ax.grid()\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# waveform dictionary structure (for step interpolated quantity)\n", "# e.g. this could be a proportional PID gain over time\n", "waveform = {\n", " \"times\": np.array([0.0, 0.15, 0.3, 0.5]),\n", " \"vals\": np.array([100, 150, 75, 100]),\n", "}\n", "\n", "# plot what it looks like \n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "ax.step(waveform['times'], waveform['vals'], where='post',\n", " color='k', linewidth=1, linestyle=\"--\", marker=\"x\", markersize=7,\n", " label=\"step interpolated waveform\")\n", "\n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Scalar quantity [units]\")\n", "ax.grid()\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we define each of the controller settings, let us define some global times, i.e. the start and end of the simulation say. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# choose some min amd max simulation times\n", "tmin = 0.0\n", "tmax = 0.25" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we define the desired control settings for the **plasma category**.\n", "\n", "These waveforms govern how the plasma current, $I_p$, is controlled—either via feedback (FB), feedforward (FF), or a combination of both.\n", "\n", "The following quantities must be specified:\n", "\n", "- **`ip_ref`**: Plasma current feedback reference waveform \\[A\\] (linear).\n", "- **`vloop_ff`**: Loop voltage feedforward reference waveform \\[V·s⁻¹\\] (linear).\n", "- **`blend`**: Control blending parameter [dimensionless] (linear):\n", " - `1` → purely FB control \n", " - `0` → purely FF control \n", " - values in `(0, 1)` → mixed FB/FF control.\n", "- **`k_prop`**: Proportional gain waveform for the FB PID controller \\[s⁻¹\\] (step).\n", "- **`k_int`**: Integral gain waveform for the FB PID controller \\[s⁻²\\] (step).\n", "- **`k_deriv`**: Derivative gain waveform for the FB PID controller (no units) (step). \n", "- **`M_solenoid`**: Mutual inductance between the plasma and solenoid (required for FF control) \\[V·s·A⁻¹\\] (step).\n", "\n", "In this particular example, we will tell the PCS to hold the plasma current constant at it's starting value. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PLASMA CATEGORY DATA\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "plasma_data = {}\n", "\n", "# plasma current feedback (FB) reference waveform \n", "# here we hold the value for the initial equilibrium constant\n", "plasma_data[\"ip_ref\"] = {\n", " 'times': np.array([tmin, tmax]),\n", " 'vals': np.array([eq._profiles.Ip, eq._profiles.Ip])\n", " }\n", "\n", "# loop voltage feedforward (FF) reference waveform (can specify the desired loop voltage on the plasma)\n", "# (here we don't use a loop voltage as we will use the FB reference waveform above)\n", "plasma_data[\"vloop_ff\"] = {\n", " 'times': np.array([tmin, tmax]),\n", " 'vals': np.array([0, 0])\n", " }\n", "\n", "# blending waveform: tells controller to use purely FB control (1), purely FF control (0), or a mixture (between 0 and 1)\n", "# we are doing pure FB control here\n", "plasma_data[\"ip_blend\"] = {\n", " 'times': np.array([tmin, tmax]),\n", " 'vals': np.array([1, 1])\n", " }\n", "\n", "# proportional gain for the FB PID controller\n", "plasma_data[\"k_prop\"] = {\n", " 'times': np.array([tmin]),\n", " 'vals': np.array([-5*4])\n", " }\n", "\n", "# integral gain for the FB PID controller\n", "plasma_data[\"k_int\"] = {\n", " 'times': np.array([tmin]),\n", " 'vals': np.array([-50*2])\n", " }\n", "\n", "# derivative gain for the FB PID controller\n", "plasma_data[\"k_deriv\"] = {\n", " 'times': np.array([tmin]),\n", " 'vals': np.array([0])\n", " }\n", "\n", "# mutual inductance between plasma and solenoid: required for FF control\n", "# not used here as we do not set a loop voltage\n", "plasma_data[\"M_solenoid\"] = {\n", " \"times\": np.array([0]), \n", " \"vals\": np.array([1.0])\n", " }\n", "\n", "plasma_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the desired control settings for the **shape category**.\n", "\n", "First, we define a function that extracts the relevant shape parameters from the equilibrium object (here referred to as `plasma_descriptors`). In this example, we extract the following quantities:\n", "\n", "- **`Rin`**: Inboard midplane radius.\n", "- **`Rout`**: Outboard midplane radius.\n", "- **`Rx`**: Radial position of the lower X-point.\n", "- **`Zx`**: Vertical position of the lower X-point.\n", "\n", "These are taken directly from the equilibrium object (without measurement noise). You could instead build synthetic diagnostics (e.g. magnetic probes) to take readings around the equilibrium (with measurement noise) and then reconstruct the plasma shape based on these. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the descriptors function (it should return an array of values and take in an eq object)\n", "# outputs (\"measurements\") from this function will be passed to the PCS later on during simulation\n", "def plasma_descriptors(eq):\n", "\n", " # inboard/outboard midplane radii\n", " RinRout = eq.innerOuterSeparatrix()\n", "\n", " # find lower X-point\n", " # define a \"box\" in which to search for the lower X-point\n", " XPT_BOX = [[0.33, -0.88], [0.95, -1.38]]\n", "\n", " # mask those points\n", " xpt_mask = (\n", " (eq.xpt[:, 0] >= XPT_BOX[0][0])\n", " & (eq.xpt[:, 0] <= XPT_BOX[1][0])\n", " & (eq.xpt[:, 1] <= XPT_BOX[0][1])\n", " & (eq.xpt[:, 1] >= XPT_BOX[1][1])\n", " )\n", " xpts = eq.xpt[xpt_mask, 0:2].squeeze()\n", " if xpts.ndim > 1 and xpts.shape[0] > 1:\n", " opt = eq.opt[0, 0:2]\n", " dists = np.linalg.norm(xpts - opt, axis=1)\n", " idx = np.argmin(dists) # index of closest point\n", " Rx, Zx = xpts[idx, :]\n", " else:\n", " Rx, Zx = xpts\n", "\n", " return np.array([RinRout[0], RinRout[1], Rx, Zx])\n", "\n", "# give these descriptors some names for use in the PCS\n", "ctrl_targets = ['Rin','Rout', 'Rx', 'Zx']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we specify the waveforms used to control these shape parameters throughout the simulation. \n", "\n", "Not all shape parameters need to be defined in units of metres—these are simply illustrative examples.\n", "\n", "For each shape parameter, the following waveforms must be specified:\n", "\n", "- **`ff`**: Feedforward (FF) waveform (linear). \n", " Not used in this example, but can be employed to impose a fixed offset on the parameter \\[m\\].\n", "\n", "- **`ref`**: Feedback (FB) reference waveform (linear). \n", " These are the desired target values of the parameter over the course of the simulation \\[m\\].\n", "\n", "- **`blend`**: Control blending waveform specifying the relative contribution of FF and FB control [dimensionless](linear).\n", " - `0` → purely FF control \n", " - `1` → purely FB control \n", " - values in `(0, 1)` → mixed FF/FB control. \n", "\n", "- **`k_prop`**: Proportional gain waveform for the FB PID controller \\[s⁻¹\\] (step).\n", "\n", "- **`k_int`**: Integral gain waveform for the FB PID controller (not used here) \\[s⁻²\\] (step).\n", "\n", "- **`k_deriv`**: Derivative gain waveform for the FB PID controller (no units) (step). \n", "\n", "- **`damping`**: Damping waveform (not used here) [dimensionless] (step). \n", "\n", "In the following, we will tell the controllers to:\n", "\n", "- hold both **`Rin`** and **`Rout`** constant. \n", "- hold **`Rx`** constant for a short period, before ramping down (to increase plasma triangularity), then ramp back up to the original value.\n", "- allow **`Zx`** to be uncontrolled." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SHAPE/DIVERTOR CATEGORY\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "shape_data = {}\n", "\n", "# waveforms for Rin\n", "shape_data[\"Rin\"] = {}\n", "shape_data[\"Rin\"][\"ff\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])} \n", "shape_data[\"Rin\"][\"ref\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.28, 0.28])} \n", "shape_data[\"Rin\"][\"blend\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([1.0, 1.0])}\n", "shape_data[\"Rin\"][\"k_prop\"] = {\"times\": np.array([tmin]), \"vals\": np.array([120])}\n", "shape_data[\"Rin\"][\"k_int\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "shape_data[\"Rin\"][\"k_deriv\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "\n", "# waveforms for Rout\n", "shape_data[\"Rout\"] = {}\n", "shape_data[\"Rout\"][\"ff\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])} \n", "shape_data[\"Rout\"][\"ref\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([1.35, 1.35])} \n", "shape_data[\"Rout\"][\"blend\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([1.0, 1.0])}\n", "shape_data[\"Rout\"][\"k_prop\"] = {\"times\": np.array([tmin]), \"vals\": np.array([170])}\n", "shape_data[\"Rout\"][\"k_int\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "shape_data[\"Rout\"][\"k_deriv\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "\n", "# waveforms for Rx\n", "shape_data[\"Rx\"] = {}\n", "shape_data[\"Rx\"][\"ff\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])} \n", "shape_data[\"Rx\"][\"ref\"] = {\"times\": np.array([tmin, 0.025, 0.075, 0.175, 0.225, tmax]), \"vals\": np.array([0.57, 0.57, 0.5, 0.5, 0.57, 0.57])} # we vary Rx\n", "shape_data[\"Rx\"][\"blend\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([1.0, 1.0])}\n", "shape_data[\"Rx\"][\"k_prop\"] = {\"times\": np.array([tmin]), \"vals\": np.array([120])}\n", "shape_data[\"Rx\"][\"k_int\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "shape_data[\"Rx\"][\"k_deriv\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "\n", "# waveforms for Zx\n", "shape_data[\"Zx\"] = {}\n", "shape_data[\"Zx\"][\"ff\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])} \n", "shape_data[\"Zx\"][\"ref\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([-1.24, -1.24])} \n", "shape_data[\"Zx\"][\"blend\"] = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])}\n", "shape_data[\"Zx\"][\"k_prop\"] = {\"times\": np.array([tmin]), \"vals\": np.array([120])}\n", "shape_data[\"Zx\"][\"k_int\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "shape_data[\"Zx\"][\"k_deriv\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.0])}\n", "\n", "shape_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the schedule of virtual circuits (VCs) in the **virtual circuits category**.\n", "\n", "For each shape parameter in the **shape category**, a VC is required to map the discrepancy between the target value (from the shape controller) and the measured value (from the equilibrium) to requests for the poloidal field (PF) coil currents. The VCs used here were obtained using the procedure described in the example notebook on VC construction. Recall that VCs are defined with units \\[m/A\\].\n", "\n", "Each VC is an array with length equal to the number of `ctrl_coils`, with entries ordered consistently with that list. When specified in the waveform dictionary, VCs are **previous-value interpolated** by the PCS class. For example, for the `Rin` controller, the VC defined at time `tmin` is used over the interval $ [t_\\mathrm{min},\\, t_\\mathrm{max})$, after which the VC defined at `tmax` is applied. Note here that we use only one VC for each shape parameter as the shape changes are not that significant. If more complex shape changes are required, then more VCs (linearised around the expected plasma shape), will be required. \n", "\n", "Note that the first entry of each shape-parameter VC is zero. This ensures that the solenoid does not contribute to shape control and is reserved exclusively for plasma current, $I_p$, control.\n", "\n", "In addition to the shape VCs, a final VC is defined for the **plasma category**. This VC maps the requested change in $I_p$ from the plasma controller to PF coil current requests (well in this case, only a request to the solenoid).\n", "\n", "Finally, optional pre-programmed feedforward (FF) requests for the PF coil currents may be specified. These allow the coil currents to be set directly if desired. In this example, no such FF requests are used, and the coil currents are driven purely via feedback control.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# VIRTUAL CIRCUITS CATEGORY (provide VCs calculated elsewhere)\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "circuits_data = {}\n", "\n", "circuits_data[\"Rin\"] = {\n", " \"times\": np.array([tmin]), \n", " \"vals\": [\n", " np.array([0, 3.1394e+04, 8.5990e+03, -8.8300e+02, -1.2460e+03, -6.6690e+03, 1.2660e+03, 2.9770e+03, 3.5350e+03, 8.3110e+03, -7.8520e+03]),\n", " ]\n", " }\n", "\n", "circuits_data[\"Rout\"] = {\n", " \"times\": np.array([tmin]), \n", " \"vals\": [\n", " np.array([0, -9.9700e+02, 1.1790e+03, 3.9700e+02, -2.0000e+02, -1.8970e+03, -5.5500e+02, -2.2640e+03, -1.5860e+03, -1.9750e+03, 5.5800e+03]),\n", " ]\n", " }\n", "\n", "circuits_data[\"Rx\"] = {\n", " \"times\": np.array([tmin]), \n", " \"vals\": [\n", " np.array([0, -3.0021e+04, 3.3980e+03, 8.6900e+03, 5.8310e+03, 1.5858e+04, 2.4000e+01, -1.2100e+02, -2.2160e+03, -9.5290e+03, 1.3570e+03]),\n", " ]\n", " }\n", "\n", "circuits_data[\"Zx\"] = {\n", " \"times\": np.array([tmin]), \n", " \"vals\": [\n", " np.zeros(11),\n", " ]\n", " }\n", "\n", "circuits_data[\"plasma\"] = {\n", " \"times\": np.array([tmin]), \n", " \"vals\": [\n", " np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),\n", " ]\n", " }\n", "\n", "# store the coil order for future reference\n", "circuits_data[\"coil_order\"] = ctrl_coils\n", "\n", "# define any feedforward coil current drives on each control coil (no drives used here)\n", "zeros_dict = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])}\n", "for coil in ctrl_coils: # linearly interpolated\n", " circuits_data[coil+\"_ref\"] = zeros_dict\n", "\n", "circuits_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cell above shows how to explicitly provide your own VCs for each shape parameter (and the plasma). \n", "\n", "The following cell shows how you can generate a VC schedule using a set of FreeGSNKE equilbiria that you might have to hand.\n", "\n", "Here, we build the VCs around the initial equilibrium generated above (it essentially builds almost the same VCs as we hard coded above). \n", "\n", "Here, the `VCGenerator` class is initialised with:\n", "- a `solver` object (i.e. the static GS solver). \n", "- a `target_calculator` to compute the plasma shape parameters (`plasma_descriptors()`)\n", "- the list of corresponding shape parameter names (`target_names`). \n", "- the full list of coils available for shape control (`coils`), only needed if `generate_fixed_schedule` is going to be used (as below)\n", "- the list of coils and targets to use in a VC computation / matrix inversion (`targets_calc`, `coils_calc`)\n", "- the list of shape targets that will use the new computed VCs, overwriting any pre-scheduled ones (`targets_ctrl`)\n", "\n", "The method `generate_fixed_schedule` then takes a list of `eq` and `profile` objects (at chosen times), computing VCs in the required waveform format above. It uses the `targets_calc`, `targets_ctrl`, `coils_calc`, and `coils` choices already given to `VCGenerator` above.\n", "\n", "The `schedule` dictionary containts the VC waveforms only for the shape parameters (and the coil ordering), these need to be combined with the plasma VC and any FF drive waveforms. \n", "\n", "Refer back to the virtual_circuits notebook for more details on VC calculations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# from freegsnke.control_loop.vc_provider import VCGenerator\n", "\n", "# # VC calculator object\n", "# VCG = VCGenerator(\n", "# solver=GSStaticSolver, # solver object\n", "# target_calculator=plasma_descriptors, # for calculating shape parameters\n", "# target_names=ctrl_targets, # names of shape parameters\n", "# coils=ctrl_coils, # coils available for shape control (only needed for generate_fixed_schedule)\n", "# coils_calc=ctrl_coils[1:], # subset of coils used for VC computation (solenoid excluded here)\n", "# targets_calc=ctrl_targets, # shapes included in the VC calculation\n", "# targets_ctrl=[\"Rin\", \"Rout\", \"Rx\"], # shapes to control (others have VCs set to zero)\n", "# vc_update_rate=0.05, # how often (in seconds) to recompute the VCs during a simulation (if initialised in the PCS class, see comment input to PCS class below)\n", "# n_vc_workers=4, # number of works used to build the VCs\n", "# )\n", "\n", "# # list of start times, equilibria, and profiles used to generate the schedule (here we use the initial condition only, i.e one set of VCs)\n", "# vc_times = [0.0]\n", "# vc_eq_list = [eq]\n", "# vc_prof_list = [profiles]\n", "\n", "# # generate schedule (shape-target VCs only -- the plasma VC and feedforward\n", "# # coil drives are not produced here and are already defined separately above)\n", "# circuits_data = VCG.generate_fixed_schedule(\n", "# times=vc_times,\n", "# eq_list=vc_eq_list,\n", "# profile_list=vc_prof_list,\n", "# tikhonov_lambda=1e-12*np.ones(len(ctrl_coils[1:])), # optional regularisation of the sensitivity matrix inversion (for calculating VCs)\n", "# verbose=True, # print output\n", "# )\n", "\n", "# # build the rest of the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "# circuits_data[\"plasma\"] = {\n", "# \"times\": np.array([tmin]), \n", "# \"vals\": [\n", "# np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),\n", "# ]\n", "# }\n", "\n", "# # define any feedforward coil current drives on each control coil (no drives used here)\n", "# zeros_dict = {\"times\": np.array([tmin, tmax]), \"vals\": np.array([0.0, 0.0])}\n", "# for coil in ctrl_coils: # linearly interpolated\n", "# circuits_data[coil+\"_ref\"] = zeros_dict\n", "\n", "# circuits_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the inputs for the **systems category**.\n", "\n", "This controller specifies operational limits for the coils listed in `ctrl_coils`, including:\n", "\n", "- **`min_coil_curr_lims`**: Minimum allowable coil currents \\[A\\] (step).\n", "- **`max_coil_curr_lims`**: Maximum allowable coil currents \\[A\\] (step).\n", "- **`max_coil_curr_ramp_lims`**: Maximum allowable coil current ramp rates \\[A·s⁻¹\\] (step).\n", "\n", "There is also the option to define feedforward (FF) perturbations (linear) to the control coil currents—similar to the FF drives available in the **virtual circuits category**. These are not used here.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SYSTEMS CATEGORY\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "systems_data = {}\n", "\n", "# define the min coil current limits and ramp rate limits\n", "systems_data[\"min_coil_curr_lims\"] = {\n", " 'times': [0.0],\n", " 'vals': [np.array([-10000, -7000, -9000, -9000, -9000, -9000, -9000, -9000, -9000, -10000, -10000])],\n", " }\n", "\n", "# define the max coil current limits\n", "systems_data[\"max_coil_curr_lims\"] = {\n", " 'times': [0.0],\n", " 'vals': [np.array([10000, 7000, 9000, 9000, 9000, 9000, 9000, 9000, 9000, 0, 0])],\n", " }\n", "\n", "# define the max ramp rate limits for the coils\n", "systems_data[\"max_coil_curr_ramp_lims\"] = {\n", " 'times': [0.0],\n", " 'vals': [1e12*np.ones(len(ctrl_coils))],\n", " }\n", "\n", "# define the control coil current perturbations\n", "for name in ctrl_coils: # linearly interpolated\n", " systems_data[name+\"_pert\"] = zeros_dict\n", "\n", "systems_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the inputs for the **PF category**.\n", "\n", "This category specifies the electrical properties, gains, and operational limits of the poloidal field (PF) coils. These inputs are used by the virtual PCS to convert coil current requests into physically consistent voltage commands.\n", "\n", "The following quantities must be provided:\n", "\n", "- **`R_matrix`**: Coil resistance array for the PF coils \\[Ω\\] (step). \n", " Defined as a time-dependent waveform and restricted to the coils listed in `ctrl_coils`.\n", "\n", "- **`M_FF_matrix`**: Coil mutual inductance matrix used in feedforward (FF) terms \\[Vs/A\\] (step). \n", " This accounts for inductive coupling between PF coils when computing FF voltage requests.\n", "\n", "- **`M_FB_matrix`**: Coil mutual inductance matrix used in feedback (FB) terms \\[Vs/A\\] (step).\n", " This accounts for inductive coupling between PF coils when computing FF voltage requests. Often it is identical to `M_FF_matrix`, but is provided separately for flexibility.\n", "\n", "- **`coil_gains`**: Feedback gain applied to each PF coil \\[s\\] (step). \n", " These gains scale the feedback voltage contributions on a per-coil basis.\n", "\n", "- **`coil_voltage_lims`**: Maximum allowable voltages for each PF coil \\[V\\] (step). \n", " These limits are enforced by the PCS when generating coil voltage commands.\n", "\n", "- **`coil_voltage_slew_lims`**: Maximum allowable voltage ramp rates for each PF coil \\[V·s⁻¹\\] (step). \n", " In this example, very large values are used, effectively disabling slew-rate limiting.\n", "\n", "All quantities are defined as time-dependent waveforms, even when constant, to allow future extension to time-varying PF system parameters.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PF CATEGORY\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "pf_data = {}\n", "\n", "# coil resistances\n", "pf_data[\"R_matrix\"] = {\n", " 'times': [0.0],\n", " 'vals': [tokamak.coil_resist[0:11]],\n", " }\n", "\n", "# coil mutual inductances (on feedforward terms)\n", "pf_data[\"M_FF_matrix\"] = {\n", " 'times': [0.0],\n", " 'vals': [tokamak.coil_self_ind[0:11,0:11]],\n", " }\n", "\n", "# coil mutual inductances (on feedback terms)\n", "pf_data[\"M_FB_matrix\"] = {\n", " 'times': [0.0],\n", " 'vals': [tokamak.coil_self_ind[0:11,0:11]],\n", " }\n", "\n", "\n", "# gains on the coils (for feedback term)\n", "pf_data[\"coil_gains\"] = {\n", " 'times': [0.0],\n", " 'vals': [0.015*np.ones(len(ctrl_coils))],\n", " }\n", "\n", "# limits on voltages in the coils\n", "pf_data[\"coil_voltage_lims\"] = {\n", " \"times\": np.array([0.0]),\n", " 'vals': [[2000, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750]],\n", " } \n", "\n", "# limits on voltage ramp rates in the coils (we set no limits here)\n", "pf_data[\"coil_voltage_slew_lims\"] = {\n", " 'times': [0.0],\n", " 'vals': [1e12*np.ones(len(ctrl_coils))],\n", " }\n", "\n", "pf_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the inputs for the **vertical category**.\n", "\n", "This category configures the feedback controller responsible for stabilising and controlling the vertical position of the plasma column using the vertical control coil (here the `P6` coil). \n", "\n", "The following quantities must be specified:\n", "\n", "- **`z_ref`**: Plasma vertical position reference waveform \\[m\\] (linear). \n", " This defines the desired vertical position of the plasma as a function of time.\n", "\n", "- **`k_prop`**: Proportional gain waveform for the vertical position feedback controller \\[s⁻¹\\] (step). \n", " This gain determines the strength of the restoring force in response to vertical displacement.\n", "\n", "- **`k_deriv`**: Derivative gain waveform for the vertical position feedback controller \\[s⁻²\\] (step).\n", " This term provides damping by responding to the rate of change of the vertical position.\n", "\n", "All inputs are defined as time-dependent waveforms, allowing for time-varying vertical control behaviour if required.\n", "\n", "In this example, we steadily increase the vertical position to before holding it constant thereafter. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# VERTICAL CATEGORY DATA\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "vertical_data = {}\n", "\n", "# plasma vertical position reference\n", "vertical_data[\"z_ref\"] = {\"times\": np.array([tmin, tmin+0.125, tmax]), \"vals\": np.array([0.0, 0.01, 0.01])}\n", "\n", "# proportional gain\n", "vertical_data[\"k_prop\"] = {\"times\": np.array([tmin]), \"vals\": np.array([0.025])}\n", "\n", "# derivative gain\n", "vertical_data[\"k_deriv\"] = {\"times\": np.array([tmin]), \"vals\": np.array([5e-7])}\n", "\n", "vertical_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we define the inputs for the **coil activations category**.\n", "\n", "This category specifies when each active poloidal field (PF) coil is enabled or disabled during the simulation. Coil activations are represented as time-dependent waveforms and are applied on a per-coil basis.\n", "\n", "For each coil listed in `active_coils`, the following entry is defined:\n", "\n", "- **`_activation`**: Coil activation waveform (step). \n", " A value of:\n", " - `1.0` indicates that the coil is active and available to the PCS.\n", " - `0.0` indicates that the coil is inactive and excluded from PCS.\n", "\n", "In this example, all coils are activated for the entire simulation interval $[t_\\mathrm{min},\\, t_\\mathrm{max}]$, using constant activation waveforms. This structure allows coils to be selectively enabled or disabled at different times in more advanced scenarios.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# COIL ACTIVATIONS CATEGORY\n", "\n", "# build the dictionary which will be passed to the virtual PCS (all keys below are required)\n", "coil_activation_data = {}\n", "\n", "# coil activation times\n", "default_coil_dict = {\"times\": np.array([tmin]), \"vals\": np.array([1.0])}\n", "for name in active_coils:\n", " coil_activation_data[name+\"_activation\"] = default_coil_dict\n", "\n", "coil_activation_data.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialise the virtual PCS class\n", "\n", "Now that the inputs have been prepared we can initialise the class and use in-built functions to view or plot the waveform data in each control category. This will build internal classes for each of the above listed controllers. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# FIRE UP THE PLASMA CONTROL SYSTEM\n", "from freegsnke.control_loop.pcs import PlasmaControlSystem\n", "\n", "PCS = PlasmaControlSystem(\n", " plasma_data=plasma_data,\n", " shape_data=shape_data,\n", " shape_control_mode=\"PID\",\n", " circuits_data=circuits_data,\n", " systems_data=systems_data,\n", " pf_data=pf_data,\n", " vertical_data=vertical_data,\n", " coil_activation_data=coil_activation_data,\n", " active_coils=active_coils,\n", " ctrl_coils=ctrl_coils,\n", " vertical_coils=vertical_coils,\n", " ctrl_targets=ctrl_targets,\n", " plasma_target=plasma_targets,\n", " # vc_generator=VCG, # optional: pass a VCGenerator (see commented-out cell above) to recompute VCs during the simulation below, else use pre-set schedule\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "PCS.CoilActivationController.data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# each controller should contain a copy of the waveform dictionaries internally\n", "PCS.PlasmaController.data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# it will also contain a dictionary of functions that has linearly (or previous value) interpolated the various waveforms\n", "PCS.PlasmaController.interpolants" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # it should also be possible to plot the data (uncomment following cell)\n", "# # green shading indiciates times that FB control is ON, yellow that FF is ON, a mix that both are ON, and white that there is no control\n", "# PCS.PlasmaController.plot_data(tmin=tmin, tmax=tmax)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the \"run_control\" method will be used internally later on but can also be called explicitly if required\n", "PCS.PlasmaController.run_control" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, suppose you wish to modify the shot setup data after initialising the PCS class. \n", "\n", "To do this, simply edit the waveforms required directly in the `.data` attribute of the controller of interest and then call `.update_interpolants()`. This has to be done to ensure that the controller refreshes any stale interpolants from a prior initialisation.\n", "\n", "Uncomment the code below to see how it works for an example in the `ShapeController`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # update the data entry with your new waveform\n", "# new_waveform = {\"times\": np.array([tmin, 0.05, 0.15, 0.35, 0.45, tmax]), \"vals\": np.array([0.28, 0.28, 0.30, 0.30, 0.28, 0.28])} # we vary Rin\n", "# PCS.ShapeController.data[\"Rin\"][\"ref\"] = new_waveform\n", "\n", "# # call the update function to refresh inteprolants (essential)\n", "# PCS.ShapeController.update_interpolants()\n", "\n", "# # plot to see new waveform\n", "# PCS.ShapeController.plot_data(targ=\"Rin\", tmin=tmin, tmax=tmax)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialise the nonlinear solver object\n", "Now choose the desired settings for the nonlinear solver object - recall the prior example notebook on this. \n", "\n", "The simulation timestep was chosen based on the vertical instability timescale in the output below (you need to always choose a value that is 5-10x smaller than this timescale) - this keeps the simulation numerically stable.\n", "\n", "This example explicitly uses coupling-based mode selection. The coupling ranking is calculated about the initial equilibrium and can evolve as the plasma moves or changes shape. The 30-mode reduction used here is therefore specific to this trajectory; for larger shape changes, validate the retained modes across representative equilibria or use `mode_selection=\"timescale\"`.\n", "\n", "The simulation timestep was chosen based on the vertical instability timescale in the output below (you need to always choose a value that is 5-10x smaller than this timescale) - this keeps the simulation numerically stable.\n", "\n", "The new **`n_linearization_workers`** argument controls the number of independent worker processes used to build both the current-response (`dIydI`) and profile-parameter-response (`dIydtheta`) Jacobian columns. `1` retains serial construction; set it higher to parallelise both the initial linearisation and later relinearisations. This setting is scoped to Jacobian construction and does not limit the inner threading of the GS solves used during ordinary evolution steps." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import nonlinear_solve\n", "\n", "stepping = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles, \n", " GSStaticSolver=GSStaticSolver,\n", " full_timestep=5e-4,\n", " plasma_resistivity=1e-7,\n", " mode_selection=\"coupling\",\n", " fix_n_vessel_modes=30,\n", " n_linearization_workers=4, # set >1 to parallelise Jacobian construction\n", " plasma_descriptor_function=plasma_descriptors,\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define FPDT simulation parameters\n", "\n", "Next set the key simulation parameters for the FPDT." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# FPDT SETUP\n", "\n", "# number of simulation time steps\n", "n = 500\n", "\n", "# simulation time step size (must be smaller than the vertical instability timescale)\n", "dt = stepping.dt_step\n", "\n", "# PCS time step (this must be whole fraction of the simulation time step)\n", "# it governs the frequency at which the PCS is called\n", "dt_PCS = dt/2\n", "\n", "# # starting time (leave as it is)\n", "# tmin = 0.0\n", "\n", "# automatically calculates time array and sets time step in solver object\n", "t_end = tmin + n*dt\n", "times = np.arange(tmin, t_end, dt)\n", "\n", "# (re-)initialise the dynamic solver with the initial eq and profiles\n", "stepping.initialize_from_ICs(eq, profiles)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before moving forward, we note that there are two time-dependent quantities that have not been set explicitly here. In this example, we hold the:\n", "\n", "- plasma resistivity constant.\n", "- plasma current density profiles constant.\n", "\n", "For modelling real plasma discharges, both of these quantities should passed as time-dependent inputs to the time-stepping loop below. For further details, see the earlier example notebooks where time-dependent resistivity and profile evolution are demonstrated.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### FPDT Simulation\n", "\n", "In the following cell, we initialise lists/arrays to store any equilbirium related data we wish to view after the simulation. See the FreeGSNKE example notebook (example03 - extracting_equilibrium_quantites) for a list of these and how to extract them. \n", "\n", "Note that some quantities can be computationally costly to extract and have not been optimised for speed yet!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# equilibrium-related data storage\n", "dynamic_psi = np.zeros((stepping.eq1.psi().shape[0], stepping.eq1.psi().shape[1], len(times))) # total poloidal flux\n", "dynamic_limiter_flag = np.zeros(len(times)) # flag if plasma is limited\n", "dynamic_psi_boundary = np.zeros(len(times)) # poloidal flux on plasma boundary\n", "dynamic_xpts = [] # list of X-point locations and associated poloidal flux\n", "dynamic_opts = [] # list of O-point locations and associated poloidal flux\n", "dynamic_currents = np.zeros((len(stepping.vessel_currents_vec),len(times))) # PF coil and vessel eigenmode currents \n", "dynamic_ip = np.zeros(len(times)) # total plasma current\n", "dynamic_shape_targets = np.zeros((len(ctrl_targets),len(times))) # shape parameters\n", "dynamic_triangularity = np.zeros(len(times)) # plasma triangularity\n", "linear_growth_rate = np.zeros(len(times)) # linearised (deformable) growth rate estimate\n", "\n", "# PCS-related data storage\n", "V_approved = np.zeros((len(ctrl_coils)+1,len(times))) # final PF coil voltages from the PCS class (passed into solver)\n", "ip_hist = np.zeros(len(times)) # integral term feeding into Plasma Category PID FB controller (at prior time step)\n", "ip_err = np.zeros(len(times)) # (error) for derivative term feeding into Plasma Category PID FB controller (at prior time step)\n", "T_err = np.zeros((len(ctrl_targets),len(times))) # proportional term feeding into Shape Category PID FB controller (at prior time step)\n", "T_hist = np.zeros((len(ctrl_targets),len(times))) # integral term feeding into Shape Category PID FB controller (at prior time step)\n", "I_approved = np.zeros((len(ctrl_coils),len(times))) # approved PF coil currents from prior time step feeding into System Category\n", "coil_resists = np.zeros((len(active_coils),len(times))) # PF coil resistances (to tell solver if a coil is switched on or off)\n", "z_current = np.zeros(len(times)) # vertical position of the plasma (average jtor position)\n", "dynamic_timings = np.zeros(len(times)) # solver runtime at each time step\n", "dynamic_jtor_norm = np.zeros(len(times)) # norm change in jtor between time steps (used to trigger relinearisation)\n", "threshold = 0.05 # relative jtor threshold (since last linearisation) above which relinearisation is triggered (here not used)\n", "\n", "# extract any initial values from the initial equilibrium\n", "dynamic_psi[:,:,0] = stepping.eq1.psi()\n", "dynamic_limiter_flag[0] = stepping.eq1._profiles.flag_limiter\n", "dynamic_psi_boundary[0] = stepping.eq1._profiles.psi_bndry\n", "dynamic_xpts.append(stepping.eq1.xpt)\n", "dynamic_opts.append(stepping.eq1.opt)\n", "dynamic_currents[:,0] = stepping.vessel_currents_vec.copy()\n", "dynamic_ip[0] = stepping.profiles1.Ip\n", "dynamic_shape_targets[:,0] = plasma_descriptors(eq=stepping.eq1)\n", "dynamic_triangularity[0] = eq.triangularity()\n", "linear_growth_rate[0] = np.max(np.real(stepping.linearised_sol.all_timescales))\n", "z_current[0] = stepping.eq1.Zcurrent()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before running the following cell (which takes a few minutes), be sure to familiarise yourself with all of the steps.\n", "\n", "During each time step:\n", "\n", "1. The primary method in the `PCS` class, `calculate_ctrl_voltages`, is called. Using “measurements” from the current equilibrium (plasma current, coil currents, and shape parameters), it computes the voltages to be applied to the evolutive solver, as well as the coil resistances (to determine whether any coils are switched off). These voltages enact the plasma control.\n", "\n", "2. The plasma current density profile parameters are placed into a dictionary for use by the evolutive solver. \n", " *Note:* they are constant here, but they can be made time-dependent at this stage.\n", "\n", "3. The evolutive solver is then called using these parameters (along with additional ones). At this stage, several choices can be made:\n", " - Specify a time-dependent plasma resistivity (usually required for resimulation of prior discharges).\n", " - Select, via `linear_only`, either a linear or fully nonlinear solution of the circuit, plasma, and Grad–Shafranov equations.\n", " - Set a relinearisation threshold when using linear mode (see example notebook 05c).\n", " - Choose not to solve the Grad–Shafranov equation fully (when using linear mode), but instead solve only for the shape parameters specified in `plasma_descriptors` (see example notebooks 05b and 05c).\n", "\n", "4. The required data are stored following successful completion of the time step.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# RUN FPDT SIMULATION\n", "for i, t in enumerate(times[0:-1]):\n", " print(\"-----\")\n", " print(f\"t = {np.round(t,5)}s (step {i+1}/{n-1})\")\n", " \n", " # start timer\n", " start_time = time.time()\n", "\n", " # initialise any historical quantities for PCS PID controllers\n", " if i == 0:\n", " ip_hist_prev = 0.0\n", " ip_err_prev = 0.0\n", " T_err_prev = np.zeros(len(ctrl_targets))\n", " T_hist_prev = np.zeros(len(ctrl_targets))\n", " I_approved_prev = eq.tokamak.getCurrentsVec()[0:11] # must use starting currents here\n", " V_approved_prev = np.zeros(len(ctrl_coils))\n", " zipv_meas = 0.0\n", " else:\n", " ip_hist_prev=ip_hist[i-1].copy() # integral term from FB PID Plasma controller at prior step\n", " ip_err_prev=ip_err[i-1].copy() # error term from FB PID Plasma controller at prior step\n", " T_err_prev=T_err[:,i-1].copy() # proportional term from FB PID Shape controller at prior step\n", " T_hist_prev=T_hist[:,i-1].copy() # integral term from FB PID Shape controller at prior step\n", " I_approved_prev=I_approved[:,i-1].copy() # approved PF coil currents from prior step (ctrl_coils only)\n", " V_approved_prev=V_approved[0:-1,i-1].copy() # approved PF coil voltages from prior step (ctrl_coils only)\n", " zipv_meas = ((z_current[i]-z_current[i-1])/dt)*dynamic_ip[i] # rate of change of vertical plasma position \n", "\n", " # # if you wish to calculate new VCs on the fly uncomment below\n", " # vcg_inputs = VCG.get_inputs_from_eq(eq=stepping.eq1, profiles=stepping.profiles1)\n", " \n", " # call the PCS class to attain PF coil voltages (on ctrl_coils and vertical coil)\n", " V_approved[:,i], ip_hist[i], ip_err[i], T_err[:,i], T_hist[:,i], I_approved[:,i], coil_resists[:,i] = PCS.calculate_ctrl_voltages(\n", " t=t, # current simulation time\n", " dt=dt_PCS, # PCS time step\n", " dt_simulator=dt, # simulator (solver) time step\n", " ip_meas=dynamic_ip[i], # measured plasma current\n", " ip_hist_prev=ip_hist_prev, # plasma controller integral term history\n", " ip_err_prev=ip_err_prev, # plasma controller error term history (for derivative term)\n", " T_meas=dynamic_shape_targets[:,i].copy(), # measured shape parameters\n", " T_err_prev=T_err_prev, # shape controller proportional term history\n", " T_hist_prev=T_hist_prev, # shape controller integral term history\n", " I_approved_prev=I_approved_prev, # approved PF currents from prior timestep\n", " I_meas=dynamic_currents[0:11,i].copy(), # measured PF coil currents\n", " V_approved_prev=V_approved_prev, # approved PF voltages from prior timestep\n", " zip_meas=z_current[i]*dynamic_ip[i], # measured vertical position x measured plasma current\n", " zipv_meas=zipv_meas, # derivative of above\n", " active_coil_resists=tokamak.coil_resist[0:12].copy(), # PF coil resistances (these are constant)\n", " verbose=False, # print some output?\n", " # vcg_inputs=vcg_inputs, # pass in eq and profiles at each time to re-calculate VCs after each time period (vc_update_rate)\n", " )\n", "\n", " # extract plasma current density profile parameters (these are constant but can be time-dependent if desired)\n", " profile_params = {\n", " \"alpha\": profiles.alpha[0:2],\n", " \"beta\": profiles.beta[0:2],\n", " }\n", "\n", " # run freegsnke over the time step with the calculated voltages\n", " stepping.nlstepper(\n", " plasma_resistivity=1e-7, # assign plasma resistivity (chosen to be constant here)\n", " active_voltage_vec=V_approved[:,i], # assign approved PF coil voltages from the PCS\n", " profiles_parameters=profile_params, # assign profile parameters\n", " custom_active_coil_resistances=coil_resists[:,i], # assign active coil resistances from PCS (tells solver if any coils are switched off)\n", " linear_only=True, # linear or nonlinear solve?\n", " target_relative_tol_currents=1e-2, # relative tolerance in the currents required for convergence\n", " target_relative_tol_GS=1e-2, # relative tolerance in the plasma flux required for convergence\n", " working_relative_tol_GS=(1e-2)/2, # tolerance used when solving GS equation, expressed in terms of the change in the plasma flux due to one timestep of evolution (must be smaller tolerance above)\n", " max_solving_iterations=20, # stop after this many iterations\n", " verbose=False, # print some output?\n", " relinearise_threshold=threshold, # if the relative jtor norm change from last linearisation is above threshold, relinearise around current equilibrium\n", " no_GS=False, # do not solve GS at each time?\n", " )\n", "\n", " # stop timer\n", " end_time = time.time()\n", "\n", " # extract and store relevant data\n", " dynamic_psi[:,:,i+1] = stepping.eq1.psi()\n", " dynamic_limiter_flag[i+1] = stepping.eq1._profiles.flag_limiter\n", " dynamic_psi_boundary[i+1] = stepping.eq1._profiles.psi_bndry\n", " dynamic_xpts.append(stepping.eq1.xpt)\n", " dynamic_opts.append(stepping.eq1.opt)\n", " dynamic_currents[:,i+1] = stepping.vessel_currents_vec.copy()\n", " dynamic_ip[i+1] = stepping.currents_vec[-1] * stepping.plasma_norm_factor\n", " dynamic_shape_targets[:,i+1] = plasma_descriptors(stepping.eq1)\n", " dynamic_timings[i] = end_time - start_time\n", " dynamic_triangularity[i] = stepping.eq1.triangularity()\n", " z_current[i+1] = stepping.eq1.Zcurrent()\n", " dynamic_jtor_norm[i+1] = stepping.relinearise_criteria\n", " stepping.linearised_sol.calculate_linear_growth_rate()\n", " linear_growth_rate[i+1] = np.max(np.real(stepping.linearised_sol.all_timescales))\n", "\n", " # print some stuff to track solve\n", " print(f\" PLASMA QUANTITIES\")\n", " print(f\" Ip = {np.round(dynamic_ip[i+1]/1000,1)} [kA]\")\n", " print(f\" Shape parameters {ctrl_targets} = {np.round(dynamic_shape_targets[:,i+1],3)}\")\n", " print(f\" Z position = {np.round(z_current[i+1],5)} [m]\")\n", " print(f\" Linearised (deformable) growth rate = {np.round(1/linear_growth_rate[i+1],3)} [1/s]\")\n", " print(f\" SOLVER STATUS\")\n", " print(f\" Relinearisation status = {np.round(dynamic_jtor_norm[i+1],3)} (threshold = {threshold}) \")\n", " \n", " # break if growth rate timescale exceeds timestep multiple\n", " if dt > linear_growth_rate[i+1]:\n", " print(\" BREAKING LOOP - INSTABILITY TIMESCALE EXCEEDS TIMESTEP SIZE\")\n", " break\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot some results\n", "\n", "Let's now plot some of the output. We can compare the evolution of the controlled parameters (plasma current, shape parameters, and vertical position) against their FB reference waveforms. In these plots:\n", "- the black dashed indicated the FB reference waveform. \n", "- the solid blue the simulated quantity.\n", "- green background shading indicates when FB control is ON. \n", "- yellow background shading indicates when FF control is ON (not used). \n", "- white background shading indicates no control is used." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PLASMA CURRENT EVOLUTION\n", "\n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "# --- references and masks ---\n", "FB_reference = PCS.PlasmaController.interpolants['ip_ref'](times)\n", "FF_reference = PCS.PlasmaController.interpolants['vloop_ff'](times)\n", "\n", "blend = PCS.PlasmaController.interpolants['ip_blend'](times)\n", "\n", "FB_mask = (blend > 0) & (np.abs(FB_reference) > 0)\n", "FF_mask = (blend < 1) & (np.abs(FF_reference) > 0)\n", "\n", "# --- shade FB regions (green) ---\n", "on_regions = np.where(np.diff(FB_mask.astype(int)) != 0)[0] + 1\n", "for seg_t, seg_state in zip(np.split(times, on_regions),\n", " np.split(FB_mask, on_regions)):\n", " if np.all(seg_state):\n", " ax.axvspan(seg_t[0], seg_t[-1],\n", " color='green', alpha=0.25,\n", " label=\"FB active\")\n", "\n", "# --- shade FF regions (yellow) ---\n", "on_regions = np.where(np.diff(FF_mask.astype(int)) != 0)[0] + 1\n", "for seg_t, seg_state in zip(np.split(times, on_regions),\n", " np.split(FF_mask, on_regions)):\n", " if np.all(seg_state):\n", " ax.axvspan(seg_t[0], seg_t[-1],\n", " color='gold', alpha=0.25,\n", " label=\"FF active\")\n", "\n", "# --- FreeGSNKE ---\n", "ax.plot(times, dynamic_ip,\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\")\n", "\n", "# --- FB reference ---\n", "ax.plot(times[FB_mask], FB_reference[FB_mask],\n", " color='k', linestyle='--',\n", " linewidth=1.5,\n", " label=\"FB reference\")\n", "\n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Plasma current [A]\")\n", "ax.grid()\n", "\n", "# deduplicate legend entries\n", "handles, labels = ax.get_legend_handles_labels()\n", "ax.legend(dict(zip(labels, handles)).values(),\n", " dict(zip(labels, handles)).keys())\n", "\n", "ax.set_ylim([7.54e5, 7.64e5])\n", "fig.suptitle(\"Plasma Current\", y=0.98)\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PLASMA VERTICAL POSITION EVOLUTION\n", "\n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "# --- FB reference ---\n", "FB_reference = PCS.VerticalController.interpolants['z_ref'](times)\n", "\n", "# --- shade FB-active region (always on here) ---\n", "ax.axvspan(times[0], times[-1],\n", " color='green', alpha=0.2,\n", " label=\"FB active\")\n", "\n", "# --- references ---\n", "ax.plot(times, FB_reference,\n", " color='k', linestyle='--',\n", " linewidth=1.5,\n", " label=\"FB reference\")\n", "\n", "# --- FreeGSNKE ---\n", "ax.plot(times, z_current,\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\")\n", "\n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Vertical position [m]\")\n", "ax.grid()\n", "ax.legend()\n", "\n", "fig.suptitle(\"Vertical Position\", y=0.98)\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SHAPE PARAMETER EVOLUTION\n", "# the ramp down in Rx here is designed to show how one might increase triangularity of the plasma (see next plot for this effect)\n", "\n", "ntarg = len(ctrl_targets)\n", "\n", "fig, axes = plt.subplots(\n", " nrows=ntarg,\n", " ncols=1,\n", " figsize=(8, 12),\n", " dpi=80,\n", " sharex=True\n", ")\n", "\n", "for j, targ in enumerate(ctrl_targets):\n", " ax = axes[j]\n", "\n", " # --- references and masks ---\n", " FF_reference = PCS.ShapeController.interpolants[targ]['ff'](times)\n", " FB_reference = PCS.ShapeController.interpolants[targ]['ref'](times)\n", "\n", " FF_mask = (\n", " (PCS.ShapeController.interpolants[targ]['blend'](times) < 1)\n", " & (np.abs(PCS.ShapeController.interpolants[targ]['ff'].derivative()(times)) > 0)\n", " )\n", "\n", " FB_mask = (\n", " (PCS.ShapeController.interpolants[targ]['blend'](times) > 0)\n", " & (np.abs(FB_reference) > 0)\n", " )\n", "\n", " # --- shade FB regions (green) ---\n", " on_regions = np.where(np.diff(FB_mask.astype(int)) != 0)[0] + 1\n", " for seg_t, seg_state in zip(np.split(times, on_regions),\n", " np.split(FB_mask, on_regions)):\n", " if np.all(seg_state):\n", " ax.axvspan(seg_t[0], seg_t[-1], color='green', alpha=0.25,\n", " label=\"FB active\" if j == 0 else None)\n", "\n", " # --- shade FF regions (yellow) ---\n", " on_regions = np.where(np.diff(FF_mask.astype(int)) != 0)[0] + 1\n", " for seg_t, seg_state in zip(np.split(times, on_regions),\n", " np.split(FF_mask, on_regions)):\n", " if np.all(seg_state):\n", " ax.axvspan(seg_t[0], seg_t[-1], color='gold', alpha=0.25,\n", " label=\"FF active\" if j == 0 else None)\n", "\n", " # --- references ---\n", " ax.plot(times[FB_mask], FB_reference[FB_mask],\n", " color='k', linestyle='--', linewidth=1.5,\n", " label=\"FB reference\" if j == 0 else None)\n", "\n", " if np.any(FF_mask):\n", " ax.plot(times[FF_mask], FF_reference[FF_mask],\n", " color='r', linestyle='--', linewidth=1.5,\n", " label=\"FF reference\" if j == 0 else None)\n", "\n", " # --- FreeGSNKE ---\n", " ax.plot(times, dynamic_shape_targets[j, :],\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\" if j == 0 else None)\n", "\n", " ax.set_ylabel(f\"{ctrl_targets[j]} [m]\")\n", " ax.grid()\n", " ax.set_ylim([np.min(FB_reference)-0.02, np.max(FB_reference)+0.02])\n", "\n", "# shared x label\n", "axes[-1].set_xlabel(r\"Shot time [$s$]\")\n", "\n", "# legend once\n", "axes[0].legend(ncol=4, loc=\"upper right\")\n", "\n", "fig.suptitle(\"Shape parameters\", y=0.995, fontsize=14)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PLASMA TRIANGULARITY EVOLUTION \n", "\n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "# --- FreeGSNKE ---\n", "ax.plot(times[0:-1], dynamic_triangularity[0:-1],\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\")\n", "\n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Triangularity\")\n", "ax.grid()\n", "ax.legend()\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also show the voltages produced by the PCS and the subsequent evolution of the currents for each of the PF coils. To display the voltage/current limits, uncomment the relevant sections in the next cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PF COIL VOLTAGES AND CURRENTS\n", "\n", "ncoils = len(active_coils)\n", "\n", "fig, axes = plt.subplots(\n", " nrows=ncoils,\n", " ncols=2,\n", " figsize=(16, 45),\n", " dpi=80,\n", " sharex=True\n", ")\n", "\n", "for j, coil in enumerate(active_coils):\n", " axV = axes[j, 0] # voltage column\n", " axI = axes[j, 1] # current column\n", "\n", " # VOLTAGES (left)\n", "\n", " # # voltage limits\n", " # if coil not in [\"P6\"]:\n", " # vlim = PCS.PFController.data['coil_voltage_lims']['vals'][0][j]\n", " # axV.hlines([-vlim, vlim], times[0], times[-1],\n", " # colors='k', linestyles='--', linewidth=1.2,\n", " # label=\"Coil limits\" if j == 0 else None)\n", "\n", " axV.plot(times[0:-1], V_approved[j, 0:-1],\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\")\n", "\n", " axV.set_ylabel(f'{coil} voltage [V]')\n", " axV.grid()\n", " if j == 0:\n", " axV.legend()\n", "\n", " # CURRENTS (right)\n", "\n", " # # current limits\n", " # if coil not in [\"P6\"]:\n", " # imin = PCS.SystemsController.data['min_coil_curr_lims']['vals'][0][j]\n", " # imax = PCS.SystemsController.data['max_coil_curr_lims']['vals'][0][j]\n", " # axI.hlines([imin, imax], times[0], times[-1],\n", " # colors='k', linestyles='--', linewidth=1.2,\n", " # label=\"Coil limits\" if j == 0 else None)\n", "\n", " axI.plot(times, dynamic_currents[j, :],\n", " color='navy', linewidth=1,\n", " marker='x', markersize=0,\n", " label=\"FreeGSNKE\")\n", "\n", " axI.set_ylabel(f'{coil} current [A]')\n", " axI.grid()\n", " if j == 0:\n", " axI.legend()\n", "\n", "# x-labels only on bottom row\n", "axes[-1, 0].set_xlabel(r'Shot time [$s$]')\n", "axes[-1, 1].set_xlabel(r'Shot time [$s$]')\n", "\n", "# column titles\n", "axes[0, 0].set_title(\"PF Coil Voltages\")\n", "axes[0, 1].set_title(\"PF Coil Currents\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given that we have passive conducting structures in the machine, we can also plot the evolution of their currents. Note how they all start from zero given that we did not initialise them in this example. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PASSIVE CURRENTS EVOLUTION\n", "\n", "# enter the names of some passives here\n", "passives_to_plot = [\"vessel_1\", \"centrecolumn_1\", \"gas_baffle_upper_1\", \"p5_case_upper_1\"]\n", " \n", "# indices of passives_to_plot in tokamak.coil_names\n", "indxs = [tokamak.coil_names.index(name) for name in passives_to_plot]\n", " \n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", " \n", "# --- FreeGSNKE ---\n", "for name, idx in zip(passives_to_plot, indxs):\n", " ax.plot(\n", " times, dynamic_currents[idx, :],\n", " linewidth=1,\n", " marker='x', markersize=0,\n", " label=name,\n", " )\n", " \n", "ax.set_xlabel(r\"Shot time [$s$]\")\n", "ax.set_ylabel(\"Currents [A]\")\n", "ax.grid()\n", "ax.legend()\n", " \n", "fig.suptitle(\"Passive currents evolution\", y=0.98)\n", "plt.tight_layout()\n", "plt.show()\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we see the evolution of the vertical instability timescale (i.e. reciprocal of the growth rate). We can see that the plasma remains stable, with the timescale only changing after the internal Jacobians are re-linearised (around 0.06 seconds)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# LINEAR INSTABILITY TIMESCALE EVOLUTION (1/GROWTH_RATE)\n", "\n", "fig, ax = plt.subplots(\n", " nrows=1,\n", " ncols=1,\n", " figsize=(8, 4),\n", " dpi=80\n", ")\n", "\n", "# plot results\n", "ax.plot(times[0:len(linear_growth_rate)], linear_growth_rate, color='navy', linewidth=1, marker=None, markersize=3, label=\"FPDT\")\n", "ax.hlines(dt, xmin=times[0], xmax=times[-1], color='black', linestyle='--', linewidth=1.5, label=f\"dt\")\n", "\n", "\n", "ax.set_ylim([1e-4, 1e-1])\n", "ax.set_yscale('log')\n", "# ax.legend()\n", "ax.grid()\n", "ax.set_xlabel(r'Shot time [$s$]')\n", "ax.set_ylabel('Linearised instability timescale [s]')\n", "fig.suptitle(\"Vertical instability timescale (linearised deformable)\", y=0.98)\n", "plt.tight_layout()\n", "plt.show()\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following cell can be uncommented to make an animation of the equilibrium in the machine over time. Might take a few seconds to run and will save the output .mp4 in the `data` directory. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # SAVE AN ANIMATION OF THE SIMULATION\n", "# %matplotlib inline\n", "\n", "# import matplotlib.animation as animation\n", "\n", "# plt.rcParams['figure.dpi'] = 90\n", "\n", "# fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", "\n", "# # Plot static wall / tokamak background\n", "# eq.tokamak.plot(axis=ax1, show=False)\n", "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2)\n", "# ax1.grid(True, which='both', alpha=0.5)\n", "# ax1.set_aspect('equal')\n", "# ax1.set_xlabel(r'Major radius, $R$ $[m]$')\n", "# ax1.set_ylabel(r'Height, $Z$ $[m]$')\n", "# ax1.set_xlim(0.05, 2.15)\n", "# ax1.set_ylim(-2.25, 2.25)\n", "\n", "# # Determine contour levels from entire dynamic_psi\n", "# min_psi = np.min(dynamic_psi)\n", "# max_psi = np.max(dynamic_psi)\n", "# levels = np.linspace(min_psi, max_psi, 40)\n", "\n", "# # --- Storage for dynamic artists ---\n", "# contour_artists = []\n", "# scatter_artists = []\n", "\n", "# # --- Update function ---\n", "# def update(i):\n", " \n", "# global contour_artists, scatter_artists\n", "\n", "# # Remove previous dynamic artists\n", "# for c in contour_artists + scatter_artists:\n", "# if isinstance(c, list):\n", "# for coll in c:\n", "# coll.remove()\n", "# else:\n", "# c.remove()\n", "# contour_artists = []\n", "# scatter_artists = []\n", "\n", "# ax1.set_title(rf\"$t$ = {np.round(times[i],3)}\")\n", "\n", "# # Main psi contours\n", "# c1 = ax1.contour(eq.R, eq.Z, dynamic_psi[:,:,i], levels=levels, alpha=0.8, cmap='viridis')\n", "# contour_artists.append(c1)\n", "\n", "# # plasma boundary\n", "# c2 = ax1.contour(eq.R, eq.Z, dynamic_psi[:,:,i],\n", "# levels=[dynamic_psi_boundary[i]],\n", "# linestyles=\"-\", colors='r', linewidths=1.4)\n", "# contour_artists.append(c2)\n", "\n", "# # adds separatrix of primary X-point if plasma limited\n", "# if dynamic_limiter_flag[i]:\n", "# c3 = ax1.contour(eq.R, eq.Z, dynamic_psi[:,:,i],\n", "# levels=[dynamic_xpts[i][0,2]],\n", "# linestyles=\"--\", colors='k', linewidths=1.4)\n", "# contour_artists.append(c3)\n", "\n", "\n", "# # X-points and O-points\n", "# sc1 = ax1.scatter(dynamic_xpts[i][:,0], dynamic_xpts[i][:,1], color='r', marker='x', s=30)\n", "# sc2 = ax1.scatter(dynamic_opts[i][:,0], dynamic_opts[i][:,1], color='g', marker='2', s=30)\n", "\n", "# # indicate which is the primary X-point more clearly\n", "# sc3 = ax1.scatter(dynamic_xpts[i][0,0], dynamic_xpts[i][0,1], color='r', marker='x', s=60)\n", "\n", "# scatter_artists.extend([sc1, sc2, sc3])\n", "\n", "# return contour_artists + scatter_artists\n", "\n", "# # --- Animation ---\n", "# num_frames = len(times)\n", "# frames_to_plot = range(0, num_frames, 5)\n", "# fps = int(len(frames_to_plot)/10)\n", "# ani = animation.FuncAnimation(fig1, update, frames=frames_to_plot, interval=1, blit=True)\n", "\n", "# # save to video (much faster & smaller than GIF)\n", "# ani.save(f\"data/animated_equilibrium.mp4\", writer=\"ffmpeg\", fps=fps)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.20" } }, "nbformat": 4, "nbformat_minor": 4 }